function LE = LEs(a0, a1, a2, b0, b1, b2, d1, d2, x0, y0, z0, N_total, N_trans)
% 计算离散映射 VL-DMSIN 的李指数
%   输入参数：
%       a0, a1, a2, b0, b1, b2, d1, d2 - 系统参数
%       x0, y0, z0 - 初始状态
%       N_total - 总迭代步数
%       N_trans - 舍弃瞬态的步数
%   输出：
%       LE - 3 个李指数，按降序排列 (LE(1) 为最大李指数)
%
% 说明：本函数利用 QR 正交化法计算李指数。函数在迭代过程中，
%       对丢弃瞬态后的每一步累计 log|Rii| 的值，最后取平均得出李指数。

% 初始化状态
x = x0;
y = y0;
z = z0;

% 初始化 Q 矩阵（3×3 单位矩阵）用于 QR 正交化
Q = eye(3);

% 初始化累计对数增量向量
sum_log = zeros(1,3);

for n = 1:N_total
    % 更新状态，计算映射
    x_new = x + y;
    y_new = sin(a0 * y^2 * sin(a1 + a2*x)) + b2 * sin(z) + d1;
    z_new = b0 * y + b1 * z + d2;
    
    % 计算当前状态下的 Jacobian 矩阵
    theta = a1 + a2*x;
    sin_theta = sin(theta);
    cos_theta = cos(theta);
    arg = a0 * y^2 * sin_theta;
    cos_arg = cos(arg);
    M1 = a0 * a2 * y^2 * cos_theta * cos_arg;
    M2 = 2 * a0 * y * cos_arg * sin_theta;
    J = [1,        1,          0;
         M1,       M2,    b2 * cos(z);
         0,       b0,          b1];
    
    % 对 Q 进行更新并 QR 分解
    Q = J * Q;
    [Q, R] = qr(Q);
    
    % 舍弃瞬态后累计 log|Rii|
    if n > N_trans
        diagR = diag(R);
        sum_log = sum_log + log(abs(diagR))';
    end
    
    % 更新状态变量
    x = x_new;
    y = y_new;
    z = z_new;
end

% 有效步数
num_effective = N_total - N_trans;

% 计算平均李指数，并按降序排序（LE(1) 为最大李指数）
LE = sum_log / num_effective;
LE = sort(LE, 'descend');

end
